Comparison between advected-field and level-set methods in the study of vesicle 

dynamics 
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Phospholipidic membranes and vesicles constitute a basic element in real biological functions. 
Vesicles are viewed as a model system to mimic basic viscoelastic behaviors of some cells, like red 
blood cells. Phase field and level-set models are powerful tools to tackle dynamics of membranes and 
their coupling to the flow. These two methods are somewhat similar, but to date no bridge between 
them has been made. This is a first focus of this paper. Furthermore, a constitutive viscoelastic law 
is derived for the composite fluid: the ambient fluid and the membranes. We present two different 
approaches to deal with the membrane local incompressibility, and point out differences. Some 
numerical results following from the level-set approach are presented. 
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Keywords: 

I. INTRODUCTION 

Vesicles are closed membranes which are suspended in an aqueous solution. The membrane is made of a bilayer 
of phospholipid molecules. These molecules have a polar hydrophilic head which points towards the solvent, and 
hydrophobic tails which point towards the interior of the membrane (Figfl]) . Both at room and physiological temper- 
atures the bilayer is a two dimensional incompressible fluid. 

Human red blood cells (RBCs) are among the simplest animal cells. They are made, like vesicles, of phospholipid 
bilayer, plus a protein network (called spectrin), known also under the generic name of cytoskeleton (Figpl. The 
internal solution of RBC is made of a hemoglobin solution (a newtonian fluid) . The RBC is devoid of a nucleus and 
organella, and this is why it may be viewed as a simple cell. Its main function consists in oxygen supply to tissues. It 
is hoped that vesicles may represent a simplistic starting point to understand viscoelastic properties, dynamics and 
rheology of bio-fluids, such as blood. 

The problem of vesicles at equilibrium (equilibrium shapes) is now fairly understood jT] . The study of vesicles under 
non-equilibrium conditions has constituted a major focus in recent years [2'J28J (see recent review [29]). Vesicles under 
flow have revealed several fascinating non-equilibrium behaviors (tank-treading, tumbling, vacillating-breathing-aka 
swinging, trembling, and so on). RBCs exhibit also similar kinds of dynamics (see review [29]). This research activity 
knows nowadays an increasing interest in the fleld of physics, mechanics, applied mathematics, engineering science, 
and so on. This field has known during the past decade analytical, numerical and experimental progresses (see review 
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FIGURE f : A cartoon showing a vesicle and the molecular structure of its membrane. Picture taken from the web site of the 
Nasa Astrobiology Institute: http://astrobiology.nasa.gov/nai/ 
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FIGURE 2: A schematic view of a red blood cell membrane. Besides the bilayer of phospholipids, there is a network of proteins 
(called also cytoskeleton) and other proteins. 



|29p. Even at zero Reynolds number, where the hydrodynamics equations are linear (the Navier-Stokes equations 
reduce to the Stokes equations), the problem remains highly nontrivial due to the fluid/membrane coupling. This 
coupling triggers nonlinearity (inherent to any moving boundary problem, even if the bulk equations are linear), 
and non-locality (a disturbance of the membrane shape at some point is felt by other regions of the membranes due 
to the fact that hydrodynamics have long range effects; very much like Coulomb interactions in electrostatics). In 
mathematical terms, the velocity field in the fluid can be integrated out in favor of a closed nonlinear integro-differential 
equations. This is the so-called boundary integral formulation based on the use of the Greens function technique [30J. 
This method has been used for vesicles by several groups [^|5J|3T]. This method, despite its efficiency and precision, 
has some limitations. For example, in a situation where the internal (or external) fluid is non newtonian, and/or if 
the non zero Reynolds number limit (a situation encountered for blood flows in veins and arteries, for example) has 
to be treated, then the resulting nonlinearities of the fluid equations clearly rule out the use of the Greens function 
techniques. We have thus to resort to other methods. 

Phase-field and level-set approaches can, in principle, handle the above problems quite naturally. Phase field 
methods have known some popularity in the scientific community, and have been used in various topics where interfacial 
phenomena are present (crystal growth, fluid/fluid interfaces, and so on). More recently, a phase-field method has 
been adapted to vesicles [7ji32j. Later, other phase-field formulations have been presented for the same system [33H55] . 

In classical interfacial problems (crystal growth, vapor/liquid phase transitions, fracture...), the total length of the 
interface can grow (or shrink) without bound. In contrast, for membranes a new important ingredient comes into 
play: the membrane is a two dimensional incompressible fluid, and thus its area does not change in the course of 
time. More precisely, not only should the total area remain constant in the course of time, but also the area must be 
conserved locally. This question is not easy to deal with in practice (a reason why only global area conservation have 
been imposed in some phase-field models [33] [34] )■ Like the phase-field, the level-set approach treats the problem 
in the same spirit, but the precise formulations differ in each case, as we shall see here. How does the level-set 
formulation compare with the phase-field one is a question which has not yet been treated in the literature, and 
especially regarding the problem of vesicles and membranes. This paper addresses this question as a first focus. 



For the sake of completeness, it is worthwhile to cite other alternative methods which have been adapted for the 
study of biological membranes. These are (i) techniques based on dynamically triangulated models ^36; or particle- 
based mesoscale solvent, multi-particle collision dynamics |37], or their combination [TU], (iv) the so-called immersed 
boundary method [35^ , used to model some features of red blood cells [W , and (v) Lattice-Boltzmann methods [30] . 

In this paper we shall draw a parallel between the phase- field and level-set methods. It will be shown that, despite 
their apparent differences, the two methods have a quite number of links. The level-set that we shall present here 
will be compared to the phase-field models previously introduced in two and three dimensions [TJ [32] ■ We shall make 
a bridge between the two methods. We shall also give here a derivation of the postulated equation in Refs. [7, 32] 
regarding the treatment of the local membrane incompressibility. Finally, we shall show that the membrane force 
derived in jT] [32] can be written as a divergence of some rank two tensor, and this will allow us to write a closed 
form for the constitutive law of the fluid/membrane system. This law exhibits a viscoelastic behavior. We shall then 
briefly discuss the numerical method used to solve the level-set equations. Some numerical results will be presented. 

The scheme of the paper is as follows. In section In] we present some essential preliminaries. In section ?? we recall 



briefly the phase-field equations for vesicles (or membranes). In section IV we introduce the level-set formulation. 
Both formulations are Eulerian methods relying on an auxiliary field to capture the interface. We adopt the notation 
<j) for the phase-field and Lp for the level-set function. We then discuss the numerical scheme used for the solution 
of the level-set model, and present some numerical results. In section^ we shall make the bridge between the two 
methods (level-set and phase-field). We shall show in that section how one can write down a constitutive law for 



the composite fiuid (fiuid+membrane) . Section VI is devoted to conclusion and discussion of some future research 



directions. Some technical details are relegated into annexes. 

II. SOME PRELIMINARIES 

Both phase-field and level-set methods have the common feature that they define the interface in an implicit 
manner; that is the phase-field also adopts the level-set notion. A significant difference exists, however. In a phase- 
field approach the phase-field values in the two fluids are forced to remain parallel (or have constant values away from 
the interface) by imposing a double well potential, while in a level-set approach the (color) function is just advected 
by the flow. Both methods have therefore potential advantages and drawbacks. Enforcing the phase-field function to 
given values and controlling the interface width by appropriate potential brings a very stable interface description. 
However this might introduce volume loss and additional terms in the potential are usually added to control this 
leak. The philosophy of phase-field in this context is therefore to give an Eulerian description of the sharp interface 
by a diffuse interface of constant and controllable width. The phase-field is a function which is only used to give 
geometrical information on the interface. By contrast, in the level-set method the auxiliary field </? is merely advected 
by the velocity field of the continuous medium, starting initially from a signed distance to the interface. There is 
no potential to keep the field stuck to a distance function (further details will be given below), and as the diffuse 
interface is typically considered to lie between two given values of ip, e.g. — e and -l-e, the interface width varies. The 
gain to move the level-set function by advection is that some mechanical information is recorded in Lp. For example 
|V(p| records the area change of interface. This makes it therefore possible to write elastic energy merely in terms of 
the auxiliary field and to obtain an attractive complex fiuid model of the fiuid-structure coupling. The non-constant 
width of interface has however to be taken care of. But this is made as post-processing each time the level-set function 
is needed to compute the distance to interface. 

III. PHASE-FIELD FORMULATION 

In this section let us recall the main results of the phase-field model introduced earlier[71 [3^. In that model the 
vesicle is described by an advected field </) which goes smoothly from —1 to -1-1 while crossing the membrane. Let —1 
represent the interior of the vesicle and -1-1 its exterior, the membrane being localized by the zero level-set of <j). In 
order to ensure a constant width for the interface, one considers the following functional 

E^ntr^ns^c[4>]= f f 1^{1 - cb'^f + ^\^ ^\A dx (1) 

when minimized in the special case of a flat interface leads to the interfacial profile given by: 

T 

(r) = tanh ( — ^ 



The above energy is called intrinsic in the sense that its role is only to define the boundary, and it should not affect 
the physics (like the forces acting on the membrane) . The physical forces are accounted for by the introduction of the 
configurational energy which is given by 

In this expression the first term stands for the curvature energy with modulus k (which is a phase-field expression of 
the Helfrich energy), while the second one is a penalization of the membrane local length variation, ^ being therefore 
a local (i.e. it depends on the given position and on time) Lagrange multiplier (we do not use C as in [7^ .32j, since 
this symbol will denote a cut-off function thereafter) . It is classical that the normal to the interface and projector on 
its tangent plane are given by 

where I is the identity tensor, while the curvature field is expressed as c = — div n. Differentiating the above energy 
with respect to the (j) gives an explicit expression of the external force 
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-K|^(c^-4G) + (FV)^c}n + ecn-hPVC ^-y^ (3) 



where G = det((ti • V)n, (^2 • V)n, n), is the Gaussian curvature and {ti,t2,n) an orthonormal trihedron. The non- 
stationnary Stokes equations are used to find the velocity field: 

EuUt - div [ri{Vu + Vu*)] - Vp = Fconfig, divu = (4) 

where e„ is a relaxation parameter (taken small enough to mimic the zero Reynolds number limit) and ry = rjouti^ + 
(j))/2 + ?7i„(l — 0)/2 the smeared viscosity (that accounts for a viscosity contrast between the interior and exterior of 
the vesicle; rjm and r]out are the viscosities of the internal and external fiuids). As the membrane is simply advected 
by flow, a simple transport equation with velocity u (i.e. d(j)/dt — where d/dt is the material derivative) should 
be used to find (j). For stability reasons, and to guarantee that the advected field function minimizes the intrinsic 
energy, and thus it represents the interface in the course of time, the phase-field is taken to obey the equation 
d4>/dt = — £,j,5E intrinsic 1^4' ; whcrc 5 / Scj) is the functional derivative and e^ a kinetic constant fixing the time scale. 
Using (IT|) one obtains the following equation: 

(j)t+u-^<j,^ e^{cjy{l - 02) + e2(A0 + c|V(/.|)) (5) 

where the term ce^jV^I has been added by "hand" in order to cancel the wall free energy of the membrane associated 
with e22|V(/)p (which is known to lead to a surface tension-like term in the force). Since the membrane does not 
exchange matter with its surrounding environment (unlike drops where molecules from bulk can migrate to surface 
and vice versa, leading to surface variation) , its surface energy is zero. The added terms guarantees the absence of a 
surface tension [7l|32] (recently another method to deal with this problem has been suggested |41|). 

Since the membrane is locally incompressible, one has to impose that the surface projected divergence of the velocity 
field must be zero. In the sharp interface picture [5, 42] the local area incompressibility can be handled by introducing 
a space and time dependent Lagrange multiplier. This amounts to writing the contribution of the membrane energy 
related to the local membrane incompressibility condition as 



E,nc= Cirm,t)dA (6) 

where the integration is performed over the vesicle membrane, and r,„ is the vector position of the membrane. The 
Lagrange multiplier C. is then determined by imposing 

(I - n (g) n) : Vm = 0. (7) 

The above expression is nothing but the divergence along the membrane (note that I — n (g) n is the projector) . The 
Lagrange multiplier ^ does not appear in the above equation which is the associated constraint, just like the pressure 
does not appear in the solenoidal constraint in incompressible hydrodynamics. However, like the pressure, we show 
in the following that ^ appears in the other equations of the model (it couples to the velocity field) . In the phase-field 



spirit [TJ [32] the idea is to define ^ everywhere (and not only along the membrane) but confine its action to the 
membrane region. We then write 

''^(r,t)|V0|dx (8) 

where dx is the volume element of the total domain. Because (j>{r) ~ tanh(r/'\/2e), it is clear that | V(?!)| is a Dirac-like 
function of width e. This implies that the energy acts in the membrane region only, as it should be. 

In |7l [32] the tension field was postulated to obey the following equation (apart from the Laplacian term which was 
included in |32| for some numerical regularization) 

^ EE ^t + u • VC = T(I - n (g) n) : Vu (9) 

at 

T is a tension-like parameter which is chosen large enough, so that we expect this to enforce (I — n ® n) : Vu ~ 0, 
that is a local membrane quasi-incompressibility. This field was built on the basis of intuition and knowledge of the 
sharp interface problem [SJIH]. It was shown in Ref. P2 by asymptotic techniques that in the limit e — >■ the above 
equation (|9l recovers the sharp interface limit (free surface divergence of the velocity field) . In (T] [32] C is said to be 
proportional to the local extension of the membrane, and indeed the right hand side of (pi is a measure of the surface 
variation, as we shall see. It will be shown in the next section how equation (|9l can be derived. 

IV. LEVEL-SET FORMULATION 

A. Introduction 

Let us now introduce the level-set approach. In the level-set formulation |43| . one still introduces a function ip 
which is negative inside the membrane and positive outside. This function is not constrained to take values between 
— 1 and +1. Rather, it is initially the signed distance to the interface, and is then advected by a transport equation 

ifit+u-Vip^O (10) 

It was shown in |44[ |45] that when u is divergence free, |V(/?| records the stretching of the interface {ip — 0}; this 
opened the way to a complete formulation in terms of if of the membrane forces. In order to localize the interface, 
one introduces a cut-off function C, i.e. a non-negative function with compact support included into [—1, 1], of unit 
mass such that -C(-) converges to the Dirac mass when e goes to zero. We used in our simulations the following 
expression C(r) = ^(1 -I- cos(7rr')) on [—1, 1], ({r) = elsewhere. This function is of unit mass. Considering ({-), this 
gives a function with support in [— e, e] which has a mass of e as readily seen by integration. This is why the scaling 
by - is necessary to maintain a unit mass and ensure the convergence of 7C(~) to the Dirac mass. By composition 
with the level-set function, we obtain -C(^) whose support is localized in the strip — e < ip < e. If i^j is a distance 
function -C( — ) converges to the Dirac mass on the curve (p ^ Q when e — > 0. In general -C(~)l^¥'l converges to this 
Dirac mass as e — >■ 0. 

Let us first consider the case where the membrane energy would depend on |V(p| only (we shall see that |V(p| is a 
measure of local surface variation). We shall call the corresponding energy elastic energy (in contrast to curvature or 
bending energy). A vesicle membrane is inextensible. In our spirit we shall allow for a finite, albeit quite small, area 
variation. If we have in mind a capsule, for example (see review |46]). then the membrane can undergo a certain local 
area variation. Therefore our energy may be used to various cases. We define the energy in the entire space provided 
that the density energy is multiplied by the appropriate localized function introduced above. We then have 

£„(^) = / E{\\/^\)-a'^)dx (11) 

Jn £ £ 



and differentiating (i.e. taking functional derivative with respect to ip) it along ( 10 1 leads to the following elastic force: 

Vp ' 



F,„= V[-B'(|V(^|)]-div 



E'i\Vp\) 



\Wp\ 



^^ ' iV^lk(^) (12) 
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Note that E stands for a constitutive law for the membrane. A trivial example is a linear law E'{r) = A(r — 1) [50J 
A is a parameter. For vesicles A is taken quite large in order to enforce quasi-inextensibility. 



A prior study [441 US] has already considered these membrane elastic forces, but not the curvature energy, which 
we consider here. We set 

£,{^) ^ f G{c{^))\Wv\-C{^) (13) 

Jn £ e 

where G is again a constitutive law for the curvature energy. A standard choice, which is compatible with [71 [35], is 
G{r) = ^r^, which is nothing but the Helfrich curvature energy density. Recall that c is the divergence of the normal 
vector, and hence c{(p) = div r^^ (note that we have the +, while in Refs [7J |31| the opposite was chosen), which 
is thus positive for convex vesicles. Our strategy to derive the force from the energy is the following: we compute 
the time derivative of the energy and equal it to minus the power of the curvature force (this is sometimes called the 
virtual power -note that a direct functional derivative can be used as well, as in Ref.[32j): 



d_ 
~dt 



£c{(p) = - Fc- udx 



From differential calculus and using the transport equation ( 10 1 solved by tp, we have 

-f.^civ) = d£c{(p){(pt) = d£c{(p){-u ■ Vip) 

where d£cif)iS) means the differential oi £c at point tp applied to the increment S. Therefore we only need to compute 
the differential of curvature energy and apply it to the increment —u ■ \/ip. We show in annex I that 



d£ciip){6) = / div 
Jn 



-G(c(^))^ + ^Pv^- (V[|V^|G'(c(^))]) 



-Ci-)Sdx 

e £ 



which by identification leads to 
Fr = div 



-G(c(^))^ + ^Pv^- (V[|V^|G'(c(^))]) 



-C(-)V^ (14) 

£ £ 



The equations for u in the level-set formulation is similar to the phase-field formulation (see (M|), with F,„ + Fc as 
source term. In fact we rather use the full Navier-Stokes equations (i.e. by including uSJu) but this induces non 
significant differences as long as the Reynolds number is small enough. Note that in this formulation there is no need 
to introduce a Lagrange multiplier and to postulate some corresponding evolution equation (as presented in the last 
section), since the stretching is encoded the (f variable. Note also that this implies (unlike in the case of the phase-field 
approach), that we have to solve a transport equation for <p, since it is under this evolution that |V(y9| records the 



membrane stretching (see next section and section V B I . Before discussing the bridges between the phase-field and 
level-set methods, we shall first present the numerical method used for the level-set method along with some numerical 
results. 

B. Numerical procedure 

The level-set model amounts to solve the following set of equations for (u, if): 

p{ut + u ■ Vu) - d\Y{-qD{u)) + Vp = F^ + Fc (15) 

divu = (16) 

ift + u-Vif^O (17) 

with appropriate initial and boundary conditions. 

As there is no additional term in the transport equation, ip does not remain a distance function for t > 0. To see 



this, we use the fact that distance function are functions with unit gradient modulus. Differentiating (10 1, one easily 
find that |V(p| verifies 

\S/p\t + u ■ V|V^| = ~\Vp\^^^ : D{u) 

therefore, starting from jVcp^j = 1 initially, we have \Vip\ = 1 for i > if and only if 



for all time. This means that the variation of u along the direction normal to the interface is zero, which is not true in 
general. As a consequence the support of the smeared delta function -C(~) will vary, making the smeared membrane 
width locally proportional to the inverse of extension. To circumvent this effect, one uses the following trick: we use 
T^-r as an approximation of the distance function to the interface, thus replacing |V(/3|-C(-) by -C,{ .J^ . ). This 
order one approximation could seem brutal, but several numerical evidences |45| proved its efficiency and precision in 
comparison with the usual renormalization trick used in level-set methods |43] . 

The numerical scheme used is a Chorin projection method on a MAC mesh, which enforces the exact divergence 
free condition, and therefore the volume conservation, at the discrete level. This is of high importance for our problem 
since change in the volume affect the shape of the niininiizer of the curvature energy. 



C. Dimensionless parameters 



The level-set main equation is given by 



p(ut + u ■ Vu) — div(r/D{u)) + Vp - 
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Let L, U, piet and T^j-ef represent characteristic length, velocity, density and viscosity scales. Accordingly we set 
X = Lx' , u — Uu' , t ~ (L/U)t' , p = Profp', V = Vrcifi P — ViciiU / L)p' , <j) — L(f>' , and e ~ Le' . Differentiating we find 



C/2 u 

ut = -j-u't,, Vu = — V'm', D{u) 



jD'iu'), diYivDiu)) = ^ diVirD'iu')), 



Vp = ryrcf-^Vy, Vip = Vip\ c{ip) = -c\p') 



In dimensionless variables (dropping the '), and for the particular case E'{r) — K{r — 1) and G{r) — ^r we get: 



Rep{ut + u ■ Vu) — div(rI?(M)) + Vp 
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?/rof A K 

are the Reynolds, Weissenberg and capillary numbers. In a shear flow one important quantity is the shear rate 7, 
from which we can express the characteristic velocity U — ^L. Substituting we finally get 

W, = -^^, C, - ^^ 

which are the dimensionless parameters from [7l[32]. Ck is a measure of bending distortion of the vesicle. The higher 
Cfc is the easiest the bending mode is (note that bending rigidity appears in the denominator. We measures the 
stretching (or dilatation) modes. Thus A has to be taken large enough in order to prevent significant membrane 
extensibility. 

There are two additional dimensionless parameters, namely the viscosity contrast A and the reduced volume (or 
reduced area in two dimensions) v, defined by 

■gin ,. Qy/irV A 

Vout' 



X 



V3D 



A3/2 



V2D 



(18) 



7r[P/27r]2 

where V (respectively A) is the vesicle volume (respectively area) and P is the perimeter in the two dimensional 
problem. In three dimensions (respectively two dimensions), v corresponds to the ratio between the actual volume 
(respectively area) over the volume (respectively area) of the a sphere (respectively circle) having the same area 
(respectively perimeter). For a sphere (or circle) v = 1, and it is less otherwise. As an example a human RBC has 
a reduced volume u ~ 0.65. The numerical results presented below correspond to two dimensions, so v will refer to 



V2D introduced in ( 18 I 



FIGURE 3: Stream lines for the tank-treading motion of a vesicle in a shear flow, at times 10 
contrast: A = 1, reduced area v — 0.7. Colors stand for iso-pressure lines. 



5, 10, 15 and 20. Viscosity 




FIGURE 4: Stream lines for the tumbling motion of vesicle in a shear flow. Viscosity contrast: 8 at times 0.2, 4, 8, 12, 16, 20, 
24, 29. Reynolds: 10""*. Colors stand for iso-pressure lines. 



D. Numerical results 



Below we present simulation results corresponding to generic situations where the vesicle makes the classical tank- 
treading or tumbling motion, and less classical motion called vacillating-breathing (an intermediate regime which has 
been described theoretically by one of us [I5j, and experimentally by Podgorsky and Steinberg). We set Re = 0.0001 
(quite close to the Stokes limit), We = 0.000025, and Ck = 0.25. Note that the fact that Ck is much larger than 
We means that at the time scale of bending modes (that time is given typically by tIj-ciL^/k) of the vesicle is larger 
than the time of the elastic mode (that time is typically given by rjrefLA. In other words the vesicle elastic response 
is quasi-instantaneous (on the scale on the physical bending response) and tries to keep the local area (perimeter in 
two dimensions) as close as possible to the initial one. Typical variations of the observed variations is of the order of 
percent. The volume is conserved with a much better accuracy (typically 1% for A^ = 64 and 0.1% for N = 128 at 
the end of the presented simulations). We first fix the viscosity contrast A = 1. We observe that the vesicle reaches a 
stationary angle (pi, as is expected. The membrane (which is fluid) undergoes a tank-treading motion. 

Figure H corresponds to the same parameters as before, but with a viscosity ratio A = 8 (the inner viscosity 8 times 
larger than in the former test), we observe a tumbling motion. It is interesting to see what happens if the external 
viscosity is lowered by a factor 8. The viscosity contrast is still the same, but in this case the external viscosity is 
lowered (rather than increasing the internal one, as done above). Tumbling motion still prevails, but here the vesicle 
is more deformed than in the previous case (the peanut shape of the vesicle is more evident (Fig|5l). This result is 




FIGURE 5: Stream lines for the Tumbling motion of vesicle in a shear flow, at time 0.2, 3, 6, 9, i2, 15, 18, 21, 24 and 27. 
Viscosity contrast: 8. Re = 8.10~^. Colors stand for iso-pressure lines. 



understood by noting that decreasing rjout is equivalent to the previous situation provided that one increase Ck (as 
well as We, while Re is reduced by a factor of 8; actually we have multiplied Re by 8 in this simulation, since we do 
not want to vary the relative effect of inertia, as we have in mind the pure Stokes limit) by the same amount (a factor 
of 4), due to the similarity properties of the two situations (actually the two situations are equivalent). Since Ck is a 
measure of bending modes, its increase allows for more flexibility of the vesicle. 

Figure [6] shows the variation of the angle of the vesicle main axis with respect to the flow direction (when the 
angle is zero this means that the vesicle is elongated and aligned along the flow direction). At low enough viscosity 
contrast the vesicle shows a tang-treading (FigJ6k) At large enough A tumbling prevails (Figlofc). We have found that 
at intermediate regime there is a vacillating-breathing regime as shown on Figj6b. Actually it can be shown from 
general arguments in two dimensions that for a fully inextensible vesicle, that a vacillating-breathing mode can not 
take place in the small deformation regime (i.e. close to a circular shape where only the second order mode is taken 
into account; the mode of a deformation about a circular shape is proportional to e"'"^ where ip is the polar angle 
and m is an integer. The second mode corresponds to m = 2). On the one hand, at larger deformation the argument 
does not hold. In addition, recent numerical studies based on the boundary integral formulation [48| did not observe 
the vacillating-breathing, but only as a transient. However, in that work the inextensibility condition was preserved 
up to a relative variation of about 10"'^ for a resolution of 128^. If the variation is larger, then the vesicle may behave 
as a capsule, and in that case it has been recently shown that the vacillating-breathing mode precedes indeed the 
tumbling one on increasing A. A systematic future study is needed in order to ascertain this more quantitatively. 

Finally we draw the phase diagram representing the three types of motion in the plane of the reduced volume 
and the viscosity contrast, as has been done in Ref.^S^. The results are reported on Fig. It] The line of transition 
towards tumbling is quantitatively close to that reported in in Ref.[S], albeit the present line is a bit lower. This 
can be attributed to conflnement (in Ref. |8] the most quantitative results are obtained from the boundary integral 
formulation in an unbounded domain). 



V. SOME BRIDGES BETWEEN THE METHODS 



A. Divergence + gradient form of the membrane forces 



Let us first show that the total membrane force (the bending and elastic contribution) can be written as divergence 
of tensor plus a gradient Let us start with the elastic force ( 12 1 and rewrite it as 



i^,„ = v[i?'(|v^|)]|v^|-C(-) 

£ e 



div 



S'dV^I) 






£ £ 



we use the following notations and elementary results from tensorial analysis: 
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0.12 




(a) Tank treading for A = 2 





(b) Vascillating breathing for A = 5.6 (c) Tumbling for A = 7. 

FIGURE 6: Angle variation depending on the viscosity ratio, for i^ = 0.8. 



(i) For two vectors a and 6, a®h is the matrix with coefficient a^bj on hne i, column j. 
(ii) This matrix is such that (a ® b)c — (b ■ c)a; thus n (g) n is the projector on n. 
(iii) For two vector fields a and b, div(a (E) b) = (Va)6 + (div b)a. 
(iv) For a scalar function / and a vector field a, V(/a) — a^ V/ + fVa. 
The B term in Fm may be transformed by taking a = jClf)^^' ^iid b — E' {\V ip\) t^^^ in (iii) to get 

B . d. (.■„v.|)X*^ic(^)) - v(lc(f )v.m|v.|,^ 



Observe that 



v(-C(-)v^) = v^ ® v(-c(^)) + -C(-)vV - v(-c(^)) ® v^ + -c(^)vV 



e e 



e e e e 



£ e ■ 
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FIGURE 7: Transition between different modes depending on tlie reduced volume for a confinement of 0.3 (confinenient is 
defined as the ratio of the vesicle effective diameter (the diameter of a circle having the same perimeter as the actual vesicle) 
over the channel width. 



(in general a®h^h®a but equality occurs when a and h are coUinear). Moreover from (ii) one has 

(V(ic(^))0V^)^ = |V^|V(ic(7)) 
e e |V(^| £ e 



and 



Thus B may be written as 



Y72 ^"^ 



V|V(^|. 



5 ^ div ( i^'(iv^l)^^ Jc(f )) - i^'(|v^|)v(|v^| lc(f )) 



and combined with A so that Fm reads 



F^ - V {i.'(|v^|)|v^|ic(f )} - div (i.'(|v^|)^^|^ic(^) 



|V^| 



(19) 



This form of F,„ is interesting since the gradient term is useless when plugged as a source term to incompressible 
(Navier-)Stokes equation (it is absorbed by the pressure gradient). Another way consists in including the gradient 
term as the divergence of a diagonal tensor: 



F„ 



div(i.'(|V,|)|V,|(I-X^)lc(f) 



(20) 



The remaining divergence term may be grouped with the divergence term of the fluid equation, turning the fluid- 
structure problem into the study of a fluid containing a viscoelastic band having the following constitutive law: 

This form is interesting since in the particular case E'{r) — Xr (which corresponds to an elastic interface with zero area 
surface at rest) we recover the stress tensor of a Korteweg fluid (setting $ = Z{ip/e) with Z' = C)- This remark has 
been used in [49J to prove existence results for fluid-structure coupling problems formulated in Eulerian coordinates. 
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Of course this remark holds also for the phase-field formulation ^T . In this case the membrane force is given by 



V(^|V0|)-divU 



\v^\ 



(21) 



Note that I — Y^^y is the projector on the tangent subspace to the vesicle membrane. 



B. Link bet^veen ^ and \Vip\ 

The first apparent difference between the level-set and the phase-field methods is that in one case (the level-set) we 
have a pure advection, while in the second one we enforce the values of </) to be ±1 inside and outside the vesicle by 
introducing a double well potential. Remember, however (see section IVB I, that tp if it is initially a distance function, 
this properties is not preserved under dynamics. One has thus usually to resort to a renormalization such as to keep 
it a distance function |43] (note that in section IVB we propose an alternative as compared to the usual method). 
Thus this procedure is somewhat similar to forcing (p in the phase-field to remain constant in the bulk phases. More 
precisely, first observe that introducing the function Z such that Z{r) = 2 J_ (^(s)ds — 1, the phase function Z{^) 
has a behaviour similar to (f>: it is equal to —1 when ip < — e, inside the vesicle, and to 1 when ip > e, outside the 
vesicle. 

Disregarding the gradient terms in the two expressions ( 19 ) and (211 of F^, which play no role as stated above, we 
remark that thanks to the chain rule we have 



V(p(E) V(f 1 ^.V'n 



VZ(f)0VZ(f) 

Wzm\ 



Thus in order to compare the formulations, regarding the membrane forces, we must compare ^ (see equation Isl) from 
the phase-field method with E'{\Vip\) (see equation 12) of the level-set method. It is easy to prove that if ip verifies 
a transport equation with divergence free velocity field, then |Viy9| verifies 

\Vip\t + u ■ W\Vip\ = \Vip\t ■ (t ■ Vu) 

and dividing by |V(/3| and multiplying by T to obtain 

(Tlog \Vp\)t + u ■ V(riog |Vv?|) = Ti • (i • Vu) 

i.e. Tlog |V(/3| verifies the same equation as the postulated [TJIHS] equation for ^ (equation [9| . Of course we can not 
conclude that ^ — Tlog |V(^|, since the velocity fields in the two methods should not coincide. Indeed the equation 
in the phase-field approach is not a simple transport equation as in level-set, and this fact affects the elastic forces and 
thus the dynamics. But up to these technical misfits, one may consider that ^ is more or less proportional to the log 
of the local extension of the membrane instead of being simply proportional to this stretching. Another consideration 
is that, roughly speaking, and up to the Ginzburg-Landau intrinsic energy, the phase-field method can be considered 
as the level-set method where we choose as constitutive law 



E{r) =Tr(logr-l) 



in the elastic membrane energy. 



C. Divergence + gradient form of the curvature force 



The level-set curvature force ( 14 1 (a similar form could be extracted in the phase-field formulation) can be put into 
divergence form. After a lengthy calculation (see Annex II) we can write F^^ under divergence form F^^ = div cr^ 
with the tensor ct^ given by 



|V(^| 



G(c(^))I 



V(/9 



® VG'(c(^)) 



-G'(c(.))(l-X^l.^- 






ic(f) (22) 



In this expression the higher order term under the divergence contains an extension of the second fundamental form 
of the surface, to the whole space. 
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D. Comparison of curvature forces with that obtained in Ref.[32] 



The above derivation was general without specifying the form of the function G. Let us take G{r) — |r^ (this is 
the Helfrich form) in ( 14 1 which now reads 



F^'^ = Kdiv 
while the curvature force of (pi) is 



c{ipf Vif 1 



2 |V(p| iVv^l 



Pv^^ (V[|V0|c(^)] 



e e 



:,AF 



K C 



2 2 



t- V(t- Vc) }V(t) 



which is equivalent to that derived in Ref.|32| (equation A13, for cq = 0, where cq is the spontaneous curvature 
considered in Ref.fSj]; recall that the definition of the curvature there has the opposite sign here). 

In two dimension. Pv(p^ (v) = {t ■ v)t thus the formula for F^^ is identical to formula F^^ of [7] IHl], up to a factor 2 
and the convention for the sign of curvature. Note that in Ref.|31| it was not realized that the force can be written as 
a divergence, as shown here. Instead the force was in the numerical study substituted by an approximate expression, 
leading to F^^ . To arrive to that expression in Refs. [3132] the property divi — was used, which is not valid in 

general. Indeed, t = -r^^ thus divt ~ — Vx0- t-^it} = —t- X .j ^ 0. In ref. [TIES] the curvature force (equation (8) 
of |31j), could be obtained from (6) only to leading order of the interface width. The present study shows in fact that 
expression (8) is in fact exact. While this remark also holds in three dimension, the algebra to prove the coincidence 
of level-set and phase field expressions is somehow cumbersome and is therefore omitted here. 



E. Comparison of curvature forces with those of |41| 



In [U] the authors found the following elastic stress for the curvature stress (we use their notation): 

T^ = -($ - divT) ®Vtp~ TD^ip 
where (considering the case Cq = 0, i.e. no spontaneous curvature) they set 



r = ac{ip) I 



|V^| 






and 



$ 



5 / x2 VtP 



^[..VV.-A.V,] 



We are going to show that up to the cut-off C term which is included in the phase-field function, this is the same as 
the transpose of a^, in the case G{r) ~ a^, which proves that the stresses are the same (recall that all geometric 
quantity express the same in terms of the phase-field or level-set equation) . As a first step we write $ in another way. 
Indeed we note that after an expansion of the curvature expression (which is c = div{\/ip/\\/ip\)) we obtain 



c(0) 



1 



|V^| 



Ai^- 



D'^ipWip ■ Wip 



One also has 

and thus 

$ = ($ 



$. 






-^c(^)^ + 2ac(0)2 = -|c(0)^ 



Wip Wip 



iv^nv^i 



Vip Vip \ 



$ 






2a 






W(p 



'^ ^ D^ipWip 



The transpose of the level-set stress is given by 
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(-D^ = a^|V^I 






+ a\I 






(Vc(</7) «) V(/3) 






(23) 



Let us show that each term in the expression of r*-^ is hidden in this expression. Consider the term TD^ip which is 
included in the last term of (23 1, if we consider only the identity / in the last projection term. The remainder is 






Now we focus on the term divT® V(^. We first compute 



(24) 



div T = div 



"^(-)(^-^^^) 






Vip Vtp \ D'^tp \7tp 



|V(p| |v^i; iv^iiv^i 



— ac{(p)' 



Vip 



Thus taking the tensorial product gives 

Vtp Vip ' 



divT«) Vv? = a I- 



|V^| |V^| 



Vc{p) (E)Vip- acip) [ I - ^^ (g) jj^ 1 D^p i ^"P ^^ ^"P 



|V(p| |V^| 



|V^| |V^| 






The first term is exactly the second term of ( |23| . The second term has the wrong sign to match ( 24 1 . But combined 
with the new expression for —^®Vp, namely 



M^(i^1^^1^\d^ 



this gives the remaining terms, up to the extra term in the LS stress 



z^VVyj (g) Vv? 



a- 



|v^|: 



which is a spherical tensor, thus not modifying the dynamics in an incompressible flow. Therefore the two curvature 
stress tensors are identical. 



VI. CONCLUSION AND FORTHCOMING WORKS 



In this paper we investigated links between the phase-field and level-set modeling of immersed elastic membranes 
subject to curvature energy. We proved that these two formulations are equivalent from a theoretical point of 
view and validated the level-set formulation by recovering known behavior of phospholipid vesicles under shear flow. 
Undergoing works concern the generalization of our model to full membrane energy for the elastic membrane. Indeed 
while phospholipid vesicles only react to local change of area and curvature, red blood cells also react to shear in the 
tangent plane to their membrane. This is due to the spectrin network underneath. Therefore there is a need for the 
modeling of the full membrane energy of an immersed interface. This question is currently under investigation. 
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VII. ANNEXES 
A. Annex I : differential of curvature energy 

Differentiating this energy with respect to ip gives 



as.ms) ^ /^G'(c(^))div (^ - ^) |v,|lc(f )^. 



G(cM)^^Jc(7) + G(cM)|V^|-^C'(7)Mx (25) 



Integration by part of the second term yields: 



G(c(^))c(^)-C(^)<5 + G(cM)^4c'(-)V'^^ + VG(c(^))-^-C(-)'5. 



The second term cancels with the last one of equation (25 1. Denoting by Pvip^ the linear projection operator on 



iv^nv^r 



we thus have 



df,M(<5) = / G'(c(^)) div (^'^"^j^^) IV^I^C(-) - G(cM)c(^)^C(-)'5 - VG(c(^)) • ^X^(^)^dx 



in \ |vifi / £ 

which thanks to the expression for c{ip) also reads 



£ £ 



|V.^| £ V- 



dfeM('5)= / G'(c((^))div 



j_^j|V^|ic(^)-div(G(cM)^)ic(-)'5d^ 



Since Pvm^ l^*^) ' ^V — the first term integrates by part into 



V[|V^|G'(c(v.))] •Pv^4V^)-^ic(-; 



|V.^| £ £• 



V^^[|Vv.|VG'(cM)].^-C(^)da; 

O |V(y3|£ £ 



where we used the symmetry property of the projection on Viy9^. Integrating by parts once again, there holds 



dE,(ip){S) = I div 
In 



-Gic{cp)) 



1 



V</3 



Vip- 



(V[|V^|G'(c(^))]) 



-C(-)<5da; 

£ £ 



B. Annex II : divergence form of curvature energy 



To get that divergence form we start from ( 14 ) that we recall for convenience 

1 



Fr = div 



-Gici^)) 






and first use the tensor identity (iii) in VA which gives 
\'ip(g>V(p 1 



Fr = div 



-G(cM)- 



\^^ (V[|V^|G'(c(^))]) 



1..^^ 



£ £ 



|V^| iv^i 



V^®Pv^^(V[|V^|G'(c(^))]) -C(t) 



£ £ ■ 

V ( ^C(f )V^) X (^-G(c(^))^ + ^^^^" (V|Vvp|G'(cM)) 
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that we denote by div A — B. Not e that x denotes a matrix-vector product in the above formula. Working on the 
second term and using (iv) of IvA) V {\C,{^)Vip) =Vip® ^^Cij) + K(f )^V = ^C'(7)V(^ «) V^ + jC(f )^V = 



ViC(f ) (8) Vv3 + lC{^)D^ip. Moreover, recall that 



|V</p| 



V|V<^|and (V(^C(p)®V(^)Pv^4u) = 



for any u using (ii) from VA Thus 



B = -G(cM)|V^|ViC(r) - G(c(^))7C(-)V|V^| + ^Pv^- (V|V^|G'(cM)) Jc(f ) 



i«?^ 



£ £ 



-G{c{^j)V (^|V(p|ic(f )) + I^Pv^- (V|V^|G'(cM)) ki 



KrV. 



From the definition of projector, 



\^±{u)= [1- 






and by using u (E) Av — {u^ v)A'^ , 



V 



|V^| 






Let us try to write i? as a divergence term minus a remainder term. To start with, 

B = -V (^G(cM)lV^I Jc(f )) + G'(cM)VcM|V^|^C(f ) + V (^^) V (|Vv.|G'(c(^))) ^Clf ) 
and computing Vc((^) leads to 



VcM=V(div^)=div 






Therefore 



S = -V(...) + |G'(cM)|V^|divv(^^)'' + v(^^)'^V(|V^|G'(cM))Uc(f) 



-V (...)+ div G'(c((p))|V(^|V 



V^ 



l«f, 






-V (...) + div ( G'(cM)|V^|V (^^ j"^ Jc(f ) J - G'(cM)|V^|V (^)^ V (ic(f ) 



but the last term is zero thanks to the expression of V ( t^^ ) with the projector on Vip-^ and the fact that V (^Clf)) 

is colinear to V(p. Writing the gradient term as the divergence of a diagonal tensor we get a first expression in divergence 
form: 



Fr = div 



-G(cM)|V^| ( I - ^^S^l + S ® Pv^-V (|V^|G'(cM)) 



|V^|^ 



|V^| 



-G'(cM)|V¥.|V 






£ £ 



We may arrange terms further by using u ® Av — [u® v)A"^ , 

\'ip(8)\'ip 



^^ -Pv.^(-) ^^ 



\V^\ 



\Vip\ 



|V^|^ 



«, = ,^«.Vi-''^^^^ 



IVv^l 



|V^| = 
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which with v — V\V f\G' {c{ip)) gives 



^.P...(V,V.,G'(c(.))). (^.V(,V.,G'(c(.)))) (l-^^ 



^' w^» l^ ^ ^ J r iv^ J + ^^ ^ ^^'^^^^» I' iv,i 



and 



T 



Finally we get the expression Fc ~ div ct^ with 



al = |g(c(^))|V^| (^I - ^^^) + (V^ ® VG'(c(^))) {l 



2 






that we can arrange like in ( 22 1 



-^'(^(^» ( ' - i^ j ^^ [' - iv^ j I i^(p (2^) 
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